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EXECUTIVE  SUMMARY 

The  objective  of  this  work  is  to  develop  a  robust  design  methodology  for  optimizing  microstructure- 
sensitive  properties  in  aircraft  components  manufactured  using  metal  forming  processes.  The  multi-scale 
forming  design  simulator  developed  as  a  part  of  this  project  provides  means  to  select  the  sequence  of  defor¬ 
mation  processes,  design  the  dies  and  preforms  for  each  process  stage  as  well  as  the  process  conditions  such 
that  a  product  is  obtained  with  desired  shape  and  microstructure.  Modeling  of  uncertainty  propagation  in 
such  multi-scale  models  of  deformation  is  extremely  complex  considering  the  nonlinear  coupled  phenomena 
that  need  to  be  accounted  for.  The  work  completed  in  the  project  addresses  key  mathematical  and  compu¬ 
tational  issues  related  to  robust  multi-scale  design  of  deformation  processes.  Our  research  accomplishments 
include  development  of  new  mathematical  models  based  on  spectral  polynomial  chaos,  support  space,  and 
entropy  maximization  techniques  for  modeling  sources  of  uncertainties  in  material  deformation  processes. 
These  models,  in  conjunction  with  multi-scale  homogenization  models,  allow  simulations  of  the  effect  of 
micros  tract  ural  variability  on  the  reliability  of  macro-scale  systems.  We  have  developed  the  first,  stochastic 
variational  multi-scale  simulator  with  an  explicit  sub-grid  model,  a  robust  deformation  process  simulator 
using  spectral  and  collocation  methods  for  simulating  uncertainties  in  metal  forming  processes.  Finally,  re¬ 
cent  developments  including  an  information  theoretic  framework  for  mndelincr  — 1 

is  summarized. 
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1  Status  of  effort 


The  most  significant  mathematical  techniques /computational  algorithms  developed  as  a  part  of  this  project 
are  explained  briefly  in  this  final  report: 

1.  Development  of  spectral  stochastic  finite  element  method  and  non- intrusive  stochastic  Gaierkin  method 
for  robust  modeling  and  design  of  deformation  processes  [1,  2t  3] 

2.  Development  of  multi-scale  sensitivity  analysis  for  designing  microstructure- sensitive  properties  in  de¬ 
formation  processes  [4,  5,  6]. 

3.  Development  of  continuum  sensitivity  analysis  technique  for  design  of  complex  three  dimensional  de¬ 
formation  processes  [7], 

4.  Development  of  stochastic  variational  multi-scale  model  with  explicit  subgrid  modeling  for  solving 
multi-scale  partial  differential  equations  (PDEs)  in  random  heterogeneous  microstructures  [8,  9]. 

5.  Development  of  maximum  entropy  techniques  for  modeling  topological  uncertainties  in  poly  crystalline 
metallic  microstructures  and  its  influence  of  homogenized  properties. [10,  11,  12] 

Detailed  information  on  mathematical  developments  are  available  from  the  above  papers  and  the  thesis  of 
students  involved  in  the  project.  Preprints  of  the  above  papers  and  thesis  copies  are  available  for  download 
at  our  web-site  http://mpdc.mae.cornell.edu 

2  Development  of  stochastic  finite  element  techniques  for  robust 
modeling  and  design  of  deformation  processes 

Two  approaches  were  developed  towards  modeling  uncertainties  in  deformation  processes.  In  the  first  tech¬ 
nique,  Spectral  Stochastic  Finite  Element  Method  (SSFEM),  a  spectral  expansion  of  the  current  configuration 
of  a  deforming  body  is  proposed  using  Legendre  chaos  expansions  to  compute  the  stochastic  deformation 
gradient  which  is  in  turn  used  to  compute  statistics  of  several  critical  quantities  in  large  deformation  analysis. 
A  stochastic  large  deformation  analysis  following  this  approach  is  presented  in  out  work  in  [lj.  The  second 
algorithm  is  based  on  a  finite  element  representation  of  the  support  space  of  the  random  variables.  This 
method  is  particularly  useful  for  capturing  instabilities  and  bifurcations  in  physical  phenomena  as  demon¬ 
strated  in  [2j.  The  support  space  representation  can  lead  to  a  non-intrusive  decoupled  as  well  as  intrusive 
coupled  formulation  for  evaluating  the  stochastic  process.  The  highlight  of  the  decoupled  approach,  called 
Noil-Intrusive  Stochastic  Galerkin  (NISG)  method,  is  that  it  can  be  directly  applied  to  presently  available 
deterministic  codes  with  minimal  effort  needed  to  compute  the  complete  probability  density  function  (PDF) 
of  a  stochastic  process.  In  NISG  method,  the  stochastic  process  is  represented  over  the  support  space  using 
piecewise  continuous  orthogonal  polynomials  in  multi- dimensional  random  variables.  The  polynomials  we 
choose  are  locally  supported  element  shape  functions  used  for  representing  functions  in  the  finite  element 
method.  The  h  and  p  convergence  characteristics  of  the  discretized  stochastic  domain  are  identical  to  spatial 
finite  elements.  The  SSFEM  and  NISG  approaches  are  described  in  this  section, 

2,1  Spectral  stochastic  finite  element  techniques  for  robust  modeling  of  defer- 
mation  processes  [1] 

This  work  introduces  algorithms  for  quantifying  uncertainty  propagation  in  finite  deformation  problems. 
The  first  algorithm  is  based  on  the  the  Spectral  Stochastic  Finite  Element  Method  (SSFEM).  A  spectral 
expansion  of  the  current  configuration  of  a  deforming  body  is  proposed  using  Legendre  chaos  expansions 
to  compute  the  stochastic  deformation  gradient  which  is  in  turn  used  to  compute  the  stochastic  analogs  of 
the  various  quantities  which  appear  in  large  deformation  analysis.  The  deformation  process  is  assumed  to 
be  quasi-static  and  inertial  effects  are  ignored.  A  total  Lagrangian  FEM  formulation  is  used  to  solve  the 


direct  deformation  problem  in  which  material  occupying  an  initial  configuration  is  deformed  to  obtain  a 
configuration  Bn+i  at  time  t  =  in+i' 

Let  X  be  a  material  particle  in  Bn  and  let  x  =  x(X  ,tn+i)  be  its  location  at  time  irl+i.  The  total 
deformation  gradient  can  be  defined  as 

F(X,tn+i)  -  70x(X,tn+1)  -  —  (2.1) 

The  total  deformation  gradient  can  then  be  decomposed  into  the  elastic  and  plastic  parts  Fr  and  Fp< 
respectively-  In  the  subsequent  text,  we  adopt  the  following  notation.  All  quantities  defined  at  the  time 
t  =  tT1  will  be  denoted  by  the  subscript  n,  while  all  the  quantities  without  any  subscript  will  refer  to  time 
t  ss  tn+ 1,  The  subscript  n  +  1  will  only  be  used  explicitly  wherever  needed  to  facilitate  understanding. 

In  the  absence  of  body  forces,  the  equilibrium  equation  at  t  =  tn+ 1  can  be  expressed  in  the  initial 
configuration  Bo  as, 

V0.  P  =  0  (2.2) 

where  Vq  denotes  the  divergence  in  Bf>  and  P  is  the  first  Piola-Kirchhoff  stress  defined  in  the  initial 
configuration.  P  is  related  to  the  Cauchy  stress  T  defined  in  the  current  configuration  as  follows: 

P  =  det.FTF~T  (2,3) 


A  hyperelastic  constitutive  law  and  classical  J2  flow  theory  for  plastic  deformation  involving  state  vari¬ 
able  based  constitutive  models  are  used.  A  real-valued  space-time  stochastic  process  t//(p,f)  with  known 
probability  distribution  can  be  written  as  a  function  w(p,  £,  8)  with  p,  /  and  8  denoting  dependence  on  space, 
time  anti  the  probability  space,  respectively.  In  subsequent  sections,  any  quantity  with  a  ^-dependence  rep¬ 
resents  a  random  quantity  and  0  would  be  referred  to  as  the  random  dimension.  In  the  spectral  stochastic 
approach,  two  different  expansions,  the  Kar lumen  Loeve  Expansion  (KLE)  anti  the  Generalized  Polynomial 
Chaos  Expansion  (GPCE)  have  been  proposed  for  representing  random  quantities.  The  classical  polyno¬ 
mial  chaos  expansion  used  in  this  work  employs  Hermite  polynomials  in  multidimensional  Ar(0. 1)  random 
variables  as  a  trial  basis  for  the  probability  space  to  represent  a  stochastic  process.  Let  us  consider  a  set  of 
independent,  identically  distributed  (iid)  Gaussian  random  variables  denoted  as  {£*  W}i=o-  The  polynomial 
chaos  representation  of  a  random  process  can  be  written  as: 

00 

ui(x,t,9)  =  a0(a:1t)ro  +  E  a*,  (a:,  (0)) 

*1  =  1 

oo  i\ 

+  E  (*, t)r2«<,W. ft, («))  +  ■■■  (2-4) 

n=]  *i=i 

where  Tv(*  )  (popularly  known  as  polynomial  chaoses)  are  Hermite  polynomials  in  (£i(0)i£2(0)t  *  *  *)*  The 
above  equation  can  be  re-written  in  a  compact  form  as  follows: 


QC 

w(x.t,0)  =  E  Wi(xtt)ij{9) 

i=0 


(2.5) 


where  there  is  one-to-one  correspondence  between  the  functionals  Fp(*  )  and  ^(V?).  In  actual  practice, 
the  random  process  is  represented  by  a  truncated  chaos  expansion  with  a  limited  number  of  terms.  The 
polynomials  {<1^(19)  }JE0  f°rm  a  complete  orthogonal  basis  for  the  probability  space.  The  orthogonality 
relationship  is  defined  by 

<  >=<  ^  >  Sij  (2.C) 

where  <  >  denotes  the  inner  product  in  the  Hilbert  support  space  of  the  random  variables  In  this 

work,  since  we  only  deal  with  uniformly  distributed  input  and  material  uncertainties,  Legendre  polynomials 
are  used  in  the  PC  expansion. 


We  now  proceed  to  develop  stochastic  formulations  for  the  deterministic  large  deformation  problem. 
In  particular,  all  quantities  are  considered  random  and  expanded  in  a  polynomial  chaos  basis.  In  all  the 
equations  which  will  be  presented,  it  is  important  to  understand  that  the  stochastic  process  represents 
an  ensemble  of  deterministic  realizations  associated  with  a  PDF.  Thus  the  stochastic  governing  equations 
represent  the  governing  equations  for  each  of  the  deterministic  realizations  in  the  ensemble.  For  instance,  the 
stochastic  equilibrium  equation  represents  an  ensemble  of  equilibrium  equations  for  each  of  the  deterministic 
realizations  comprising  the  stochastic  process.  Thus  all  the  constraints  and  law's  which  are  valid  in  the 
deterministic  problem  (e,g.  objectivity  of  constitutive  law's,  detFe  >  0,  detFp  =  1)  are  equally  valid  in  the 
stochastic  problem. 

We  start  by  making  an  assumption  that  the  initial  undeformed  configuration  of  the  body  is  completely 
deterministic.  The  primary  variable  of  interest  in  stochastic  large  deformation  analysis  is  the  displacement 
which  can  be  expanded  as  follows: 

M 

un+](x,0)  =  y^Mn+i,(a;)$i(g)  (2.7) 

i=0 

where,  we  use  a  truncated  chaos  expansion  to  represent  the  random  variable. 

Computing  the  stochastic  deformation  gradient  The  stochastic  large-deformation  motion  of  a  body 
is  considered  in  this  section.  Let  the  initial  deterministic  configuration  of  the  body  be  denoted  as  By.  We 
define  a  stochastic  motion  x(X,  ,  0}  which  deforms  the  body  from  its  initial  configuration  to  a  final 
configuration  Bn+i(Q)  (Fig,  1).  The  location  of  the  material  point  xn+i(0)  in  I?ri  +  i((9)  is  thus  defined  by 

xn+1(0)  =  x(X,fn+]T0)  (2.8) 

In  the  configuration  £JTl+](f?),  the  coordinates  of  the  different  points  will  be  defined  by  the  same  random 
variables  which  define  the  motion.  The  stochastic  position  coordinates  of  the  final  configuration  can  be 
represented  as  xn+i  =  xn+it  $*(#).  Since  these  position  coordinates  completely  define  the  configuration  of 
the  body,  we  can  denote  B„+i{0)  as  Bn+it4q,  where  Bri+ 1(  denotes  a  virtual  configuration  defined  by  the 
points  This  configuration  may  not  have  any  physical  interpretation  but  is  only  defined  to  simplify 

the  understanding  of  the  stochastic  formulation. 

Using  the  above  definitions,  the  stochastic  deformation  gradient  at  time  fT1+1  can  be  expressed  as 

F(0)  =  V0x(X,  tn+l ,  8)  -  9a:(^+]’g)  (2.9) 

To  ensure  that  the  motion  is  well-defined,  it  is  required  that  detF[0)  >  0.  Thus  the  randomness  in  a 
meaningful  motion  will  propagate  in  such  a  manner  that  the  above  constraint  is  satisfied.  For  a  discretized 


and  parameterized  domain  such  as  a  finite  element  mesh,  Eq.  (2.9)  can  be  expressed  as  Follows: 


F(0)  =  P{xn+1,e)Q-\X)  (2.10) 

where, 

PaB^n+1,0)  =  (2-11) 

and 

OXa 

Q«0(X)  =  j—  (2.12) 

where  7}p  denotes  the  natural  coordinates  in  the  parameterized  domain  and  a,0  =  1 7...,nsd  where  nsd 
denotes  the  number  of  space  dimensions  of  the  problem.  Now,  we  can  expand  x*+1(0)  in  its  polynomial 
chaos  basis  as  follows 

<+lW  =<+!*<(«  (2-13) 

and  as  a  result  P  can  be  written  as 

P(Xn+ 1,9)  =  Pi(®T,+l,)$i(fl)  (2.14) 

where  each  P,  denotes  a  nsd  x  nsd  matrix  given  by 

a*  {*»+!,)  =  (2.i5) 

For  Rotational  simplicity,  the  dependencies  0  and  x  are  dropped  henceforth.  Using  Eqs.  (2.9),  {2. 12) 
and  (2.14),  the  deformation  gradient  can  now  be  written  as  follows: 

F  =  F,$i  =  {P,$1)(Q)"1  (2.16) 


Thus,  Ft  can  be  obtained  as 

Pi  =  P,Q-'  (2.17) 

Each  Fi  in  the  above  equation  can  be  interpreted  as  the  deformation  gradient  for  the  motion  which  maps 
the  initial  configuration  Bq  to  each  of  the  virtual  configurations  B%.  The  other  quantities  which  appear  in 
the  principle  of  virtual  work  statement  can  now  be  derived  from  F. 

Linearized  stochastic  virtual  work:  Having  computed  F ,  we  now  proceed  to  formulate  the  linearized 
virtual  work  equation  for  the  stochastic  finite  deformation  problem.  The  linearized  stochastic  virtual  work 
can  be  expressed  as  follows: 

I  dP.~dV0  =  -  /  P.^L  dV0+  [  t.udAo  (2.18) 

Jb„  ox  jb„  ox  j  r 

Expanding  the  random  variables  in  their  polynomial  chaos  expansions  (PCE),  we  obtain: 


clP'-^X^JdVo  =  ' 


Out  . 
P, — 

1  ax  1  J 


ti.UiQ&j  dAo 


(2.19) 


All  the  quantities  in  the  above  expression  except  BoydV^dA^  are  stochastic  in  nature.  The  stochastic 
quantities  are  computed  in  a  similar  manner  as  in  the  deterministic  problem  using  the  elementary  operations 
for  random  quantities  described  earlier.  In  the  process  of  computing  these  quantities,  analogous  expressions 
for  the  stress  and  strain  measures,  consistent  moduli,  elastic  and  plastic  deformation  gradients  as  well  as  the 
stress  update  algorithm  need  to  be  considered. 

Next  we  consider  the  term  in  Eq.  (2.19)  and  noting  that  we  use  a  total  Lagrangian  formulation,  we 
obtain  the  following: 

=  HAui  =  dFi 


where  II  is  a  suitably  defined  operator. 


(2.20) 


Fiber  angle  (radians):  0.215  0.225  0.235  0.245  0.255  0.265  0.275 
Figure  2:  One  of  the  realizations  in  the  KL  expansion  of  the  fiber  orientation  (Problem  3). 


On  substituting  Eq.  (2.20)  in  Eq.  (2.19)  and  discretizing  the  space  domain  using  finite  elements  and 
performing  the  integration  over  the  initial  configuration,  we  arrive  at  an  expression  of  the  form: 

^  Bk*k  (2.21) 

Finally  on  considering  a  Galerkin  projection,  a  linear  system  of  equations  is  obtained  for  the  unknowns  Aitj 
given  by 

K,Auj  <  >=  Bt  <  >  (2.22) 

Thus  using  the  definition  of  F  and  the  various  operations  on  random  matrices  and  scalars  a  complete 
stochastic  formulation  of  the  hyperelastic- viscoplastic  finite  deformation  problem  can  be  obtained. 

Numerical  Example:  Effect  of  uncertain  fiber  orientation  on  the  response  of  a  composite 
material 

In  this  problem,  we  examine  the  effect  of  uncertain  fiber  orientation  on  the  response  of  a  typical  aerospace 
component  made  of  a  composite  material.  The  2D  section  of  a  nozzle  flap  is  shown  in  Figure  2.  Plane  strain 
conditions  are  assumed.  The  total  length  is  taken  as  1  mm.  Other  dimensions  are  assumed  accordingly.  The 
flap  is  subjected  to  an  aerodynamic  force  (uniform  pressure)  on  its  back.  The  pin-eyes  are  assumed  rigidly 
fixed.  The  material  is  assumed  to  be  orthotropic- hyperelastic.  An  orthotropic  Neo-Hookean  strain  energy 
function  is  considered  for  the  problem  which  leads  to  a  non-constant  Lagrangian  elasticity  tensor  given  by 

CK  =  XdetF(2detF  -  l)C~l  ®  C~l  4-  2[n  +  XdetF(detF  -  1)]G 

+  fry  a  ®  a  @  a  @  a  +  4  f3(a  I  +  /  &  a  ®  a)  —  a  A  (2.23) 

The  vector  a  denotes  the  unit  fiber  axis.  For  the  2D  case,  it  can  be  written  as  (cos#,  sin0,G).  For  the 
stochastic  problem,  we  consider  the  fiber  orientation  to  be  uncertain  with  the  initial  fiber  angle  defined  by 
a  KL  expansion  using  a  covariance  kernel  An  exponential  covariance  kernel  was  assumed  which  is  given  by 

=  (t02erp(-^)  (2.24) 

where  r  is  the  distance  between  the  points  Pj  and  p2.  b  is  the  correlation  length  and  is  assumed  to  be  1  mm 
and  a  =  0.3.  Thus,  the  fiber  orientation  can  be  written  as, 

00 

0(P)  =  X>\/VMP)  (2-25) 

1 

where  the  mean  fiber  orientation  is  assumed  to  be  along  the  x-axis.  We  assume  that  the  flap  is  composed  of 
short  fibers  which  justifies  the  above  model  for  considering  heterogeneity  in  the  fiber  orientation.  For  this 
problem,  we  consider  terms  up  to  i  =  2  in  the  KL  expansion.  A  comparatively  large  uniform  pressure  of 
0.002  N/mm2  is  applied  at  the  back  of  the  flap.  An  order  3  PCE  is  used  leading  to  10  terms  in  the  expansion. 
Figure  2  shows  one  of  the  realizations  of  the  heterogeneity  in  the  fiber  orientation.  The  stochastic  mean 
solution  shows  a  considerable  deviation  from  the  deterministic  solution .  The  mean  and  standard  deviation 
of  the  equivalent  stress  are  plotted  in  Fig.  3. 
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Figure  3:  Mean  (top)  and  standard  deviation  (bottom)  of  the  equivalent  stress  using  order  3  PCE  (Problem 
3). 

2.2  Non- intrusive  stochastic  Galerkin  scheme  for  robust  modeling  of  deforma¬ 

tion  processes  [2] 

In  this  method,  the  stochastic  process  is  represented  over  the  support  space  using  piecewise  continuous 
orthogonal  polynomials  in  multi- dimensional  random  variables.  The  polynomials  we  choose  are  the  locally 
supported  element  shape  functions  of  a  specific  order  used  for  representing  functions  in  the  finite  element 
method.  The  space  0  spanned  by  these  random  variables  is  discretized  using  disjoint  finite  element  subdo¬ 
mains  leading  to  a  new  discretized  space  Sh  where  h  is  an  element,  size  parameter.  This  approximation  can 
be  denoted  as 

g(x,t,6)  =  2,-)  «  gh(x,t,  -4n)  =  gh(x 

Thus  the  stochastic  process  can  be  represented  using  the  basis  functions  as 

nodes 

gh(x*t,£)=  Y,  «?*(«)  (2-26) 

£=1 

where  are  the  locally  supported  basis  functions  and  gt  are  the  corresponding  nodal  values.  A  scheme  For 
discretization  of  a  one  dimensional  (AT  =  1)  pdf  is  shown  in  Fig. 4 

It  was  shown  in  [13]  that  for  functions  involving  high  non  linearities  and  critical  points  the  GPCE  fails 
to  perform  satisfactorily.  On  the  other  hand,  a  piecewrise  representation  is  ideally  suited  for  representing 
discontinuities  and  nonlinearities  in  the  processes.  As  mentioned  before,  finite  deformation  processes  involve 
considerable  degrees  on  nonlinearities  at  every  stage.  Thus  we  resort  to  the  third  approach  -  a  piecewise  rep¬ 
resentation  using  finite  element  basis  functions  for  representation  of  the  random  large  deformation  processes 
in  this  work.  This  is  also  advantageous  from  a  practical  point  of  view  since  in  converting  a  deterministic 
finite  element  based  deterministic  large  deformation  code  to  simulate  stochastic  deformation  processes,  these 
basis  functions  are  available  from  other  parts  of  the  code  and  can  be  readily  used  for  representing  random 
functions. 

The  piecewise  representation  of  the  stochastic  process  which  was  discussed  in  the  previous  section,  can 
lead  to  a  non-intrusive  decoupled  as  well  as  intrusive  coupled  formulation  for  evaluating  the  stochastic 
process.  The  two  formulations  are  explained  with  an  example  of  a  solution  of  a  linear  system.  We  consider 
the  evaluation  of  a  random  vector  u(sc,£)  by  the  solution  of  a  linear  system  which  is  defined  by  a  random 
matrix  A(£)  and  a  random  vector  J>(£)  as  shown  below.  (The  superscript  h  has  been  dropped  for  notationai 


FE  Grid 

Figure  4:  Finite  Element  discretization  scheme  for  a  1  dimensional  pdf 


simplicity) 

A(x,£)u(x,$)  =  b(x,£)  (2.27) 

Intrusive  Coupled  system  To  solve  the  system  using  the  coupled  formulation,  we  first  represent  the 
above  process  in  a  N  dimensional  piecewise  polynomial  basis  (equivalent  to  N  nodes  in  the  mesh) 

(2.28) 

where  iyj<  k  —  1  * , .  Ar.  Thus  the  solution  of  the  random  vector  u  involves  computing  the  coefficients  u;  in  the 
polynomial  representation  from  the  known  values  of  -4*  and  bk-  This  is  done  by  projecting  (inner  product) 
the  above  equation  into  each  orthogonal  basis  function  0/  leading  to 

A{Uj  <  >=  bk  <  >  (2.29) 

where  the  inner  product  is  defined  as  follows 

<  *(€)y($)  >=  f  zyf(£)d£  (2.30) 

Je 

Using  the  orthogonality  relation  of  the  polynomials  basis  functions  in.  Eq.  2.29  we  get 

AiUj  <  >=  bh  <  >  Su.  —  bt  <  >  (2-31) 

Thus  evaluation  of  the  Uj* s  involves  solution  of  a  coupled  system.  The  coupled  system  involves  N  x  ndof 
equations.  Moreover,  it  needs  explicit  representations  of  the  stochastic  processes  a(£)  and  fr(£)  in  terms  of 
the  basis  functions  which  may  not  be  readily  available  and  may  involve  more  intrusive  computations  and 
change  of  the  overall  structure  of  the  underlying  deterministic  code. 

Decoupled  system  *  NISG  formulation 

In  the  non- intrusive  decoupled  scheme,  the  output  stochastic  process  is  constructed  using  deterministic 
function  evaluations  at  an  optimal  number  of  points  defined  in  the  input  support  space.  Considering  the 
linear  system  (Eq,  2.27)  again  at  specific  points  £  of  the  support  space  we  get 

=  (2  32) 

As  can  be  observed,  each  function  evaluation  of  A,u  and  b  takes  place  at  fixed  points  £  in  the  random  space. 
This  can  be  equivalently  expressed  as 


i4j(®)ui(*)  =  Mi) 


(2.33) 


The  above  system  can  be  solved  for  m  at  predetermined  points  of  the  support  space  which  can  be  used  to 
construct  the  solution  u(a;,£).  The  procedure  is  explained  later  in  the  section. 

The  specific  function  evaluation  points  are  chosen  from  the  consideration  that  in  typical  problems  in¬ 
volving  random  variables,  the  quantities  of  interest  are  the  probability  density  function,  and  the  statistical 
moments  of  u.  The  pth  statistical  moment  can  be  calculated  as  follows: 


MP=  f 

J© 


(u(Bj,£))pf(Z)d£ 


(2.34) 


where  net  is  the  number  of  elements  in  the  support  space  while  O'1  represents  the  local  element  domain. 
Using  a  Gauss  quadrature  integration  scheme  with  n  integration  points  per  element,  the  above  expression 
can  be  approximated  as: 

?! 

Mp  w  Mp  =  y^y2wiM:E'fi£-))f7(0  (2.35) 

e=i t-1 

where  denote  the  abscissae  of  the  integration  points  while  Wi  denote  the  respective  weights. 

The  above  equation  can  be  equivalently  expressed  as: 

net  n 

Mp  -  (2.36) 

e=l  i=l 


where  t)  denotes  a  deterministic  space-time  dependent  process  at  each  integration  point  of  the  support 
space. 

Tims  evaluation  of  u,  is  completely  deterministic  and  can  be  achieved  using  a  readily  available  determin¬ 
istic  solution  scheme  and  hence  the  non- intrusive  terminology.  To  obtain  the  pdf  of  u  denoted  as  /„(£},  we 
first  obtain  the  values  of  u  at  the  nodal  points  of  the  discretized  random  space  from  the  integration  point 
values  U{.  This  can  be  achieved  using  a  standard  least  squares  smoothing  scheme,  /„(£)  can  then  be  ob¬ 
tained  by  generating  random  Monte  Carlo  samples  from  the  pdfs  of  the  input  random  variables  £ ,  followed  by 
evaluations  of  u  at  those  points.  It  should  be  noted  that  this  sampling  does  not  require  an  explicit  function 
evaluation  but  is  generated  from  an  interpolation  of  earlier  evaluations  using  the  appropriate  interpolation 
order. 

The  developed  scheme  for  discretization  of  the  stochastic  support  space  closely  follows  the  finite  element 
discretization  of  spatial  domains  in  various  applications  and  inherits  all  its  properties.  The  the  h  and  p 
convergence  characteristics  of  the  discretized  domain  are  identical  to  the  case  of  spatial  finite  elements. 

The  procedure  for  NISG  analysis  of  large  deformation  processes  can  be  summarized  as  follows* 


1.  Determine  the  pdf  of  the  input  uncertainties  in  terms  of  independent  random  variables  {£,  and 
compute  the  joint  pdf  /(C). 

2.  Determine  the  dimensionality  of  the  support  space.  Based  on  this  dimension  discretize  the  support 
space  using  finite  elements. 

3.  FOR  j=l  to  number  of  elements  in  the  grid 

FOR  i— 1  to  number  of  local  element  integration  points 

Compute  the  values  of  the  random  inputs  at  the  support  space  integration  point 
Solve  the  deterministic  large  deformation  problem  with  the  given  inputs  and  compute  u*  -  the 
output  of  interest. 

END 

END 

4.  Compute  u  from  the  values  of  Ui  using  least  squares  smoothing. 


f 


Figure  5:  Initial  configuration  (left),  final  configuration  with  no  initial  voids  (center)  final  configuration  with 
an  initial  void  fraction  of  0.03. 

5.  Compute  the  moments  as  required  from  the  n,  values  as  described  in  Eq.  2.35 

Thus,  the  NISG  approach  involves  deterministic  function  evaluations  at  the  integration  points  of  the 
discretized  support  space  of  random  variables  which  essentially  translates  to  a  decoupled  approach  for  eval¬ 
uating  the  additional  degrees  of  freedom  resulting  from  the  randomness  in  the  problem.  The  computational 
gains  attained  by  this  decoupling  are  significant.  Also  these  gains  can  offset  to  a  great  extent  the  issue  of 
slower  convergence  rate  compared  to  GPCE  by  the  selection  of  a  finer  grid  for  the  support  space. 

NISG  Example  -  Stochastic  estimation  of  die  underfill  caused  by  material  porosity 

This  problem  studies  the  effect,  of  a  random  voids  in  the  design  of  flashless  closed  die  forging  processes. 
The  material  is  chosen  to  be  an  Fe-2%  Si  alloy  at  1273  K.  A  hyperbolic  sine  flow  function  is  assumed.  The 
die- workpiece  friction  coefficient  was  taken  as  G.L  For  an  initial  billet  (dimensions  -  radius  .9339mm,  height 
1  mm)  with  no  voids  having  the  same  volume  as  the  die  (2.74mm3},  the  deterministic  flash  less  closed  die 
forging  process  is  shown  in  Fig.  5,  The  corresponding  final  state  for  the  deterministic  process  with  an  initial 
void  fraction  of  0.03  in  the  billet  is  also  shown  in  the  same  figure.  In  the  regions  where  there  are  compressive 
stresses,  the  void  fraction  decreases  leading  to  a  decrease  in  volume.  The  converse  happens  in  regions  of 
tensile  stresses.  As  a  result  an  underfill  is  observed  in  the  process. 

The  NISG  technique  provides  a  robust  way  to  estimate  the  statistics  of  the  extent  of  die  underfill  as  a 
result  of  a  random  distribution  of  voids  in  the  billet.  The  initial  void  distribution  is  derived  from  an  assumed 
exponential  correlation  kernel  using  a  Karhunen  Loeve  expansion  given  by  : 

rc(Pi,0,p2,0)  =  a2exp(-^)  (2.37) 

where  r  is  the  distance  between  the  points  pt  and  p2 ,  b  is  the  correlation  length  assumed  to  be  5  mm  and 
a  ”  0,1.  As  in  the  earlier  example,  a  truncated  KL  expansion  sufficient  to  represent  the  heterogeneities  for 
the  stochastic  void  fraction  can  be  written  as, 

2 

/(p)-/o(l+2>V'5T/i(p))  (2-38) 

1^1 

where  /0  —  0.Q3  is  the  mean  void  fraction,  A  9x9  grid  was  used  for  computing  the  statistics.  The  mean 
underfill  was  estimated  to  be  0,046976mm3  with  a  standard  deviation  of  0,0022mm3.  Using  a  10x10  support 
space  grid  the  mean  underfill  was  found  to  be  0,046 187mm3  with  a  standard  deviation  of  0.0023mm3.  The 
complete  pdf  of  the  underfill  observed  is  plotted  in  Fig, 6 
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Figure  6:  Probability  distribution  of  the  underfill  observed  due  to  material  porosity  in  the  flashless  closed 
die  Forging  problem, 

3  Development  of  multi-scale  sensitivity  analysis  for  designing 
microstructure-sensitive  properties  in  deformation  processes 

We  have  developed  a  versatile  multi-scale  design  simulator  for  tailoring  properties  of  materials  by  altering 
die  and  preform  shapes  in  industrial  forming  processes.  The  problem  is  a  significant  computational  challenge 
since  a  multi-scale  design  problem  for  modeling  3D  processes  approaches  a  billion  degrees  of  freedom.  We 
have  recently  acquired  a  128  Intel  processor  Linux  super  computing  cluster  with  funding  from  Department 
of  Defense  (DURIP-05)  for  solving  such  problems  in  materials-by-design,  The  mathematical  developments 
are  briefly  explained  in  this  section. 

Consider  a  macroscopic  material  point  and  let  it  be  associated  with  the  underlying  microstructure  M. 
Assume  that  the  response  of  any  crystal  of  the  polycrystal  is  determined  by  its  orientation  H *  which  is 
the  rotation  relating  the  crystal  lattice  frame,  e*,  to  a  sample  reference  frame  e*  as  e%  ~  R  e*.  The  ODF 
(orientation  distribution  function),  represented  as  A(r),  describes  the  crystal  density  over  the  fundamental 
region,  A  re-orientation  map,  f,  is  associated  with  the  ODF  as,  *4(rj)  —  —  *4(s,f).  The 

representation  of  the  ODF  given  by  A(r ,  t)  is  Eulerian  and  ^4(s,  t)  is  Lagrangian. 

The  evolution  of  a  poly  crystal  generally  demands  the  numerical  solution  of  an  ODF  conservation  equation. 
The  conservation  equation  for  the  ODF  in  the  Lagrangian  description  can  be  written  as: 

I  (^A(a,t)  J(a,t)  -  ,4(s,0))  do  =  0  (3.1) 

where  dv  is  the  volume  element  on  the  reference  fundamental  region,  J(s,  f)  =  del (Vf  (s,  i ))  is  the  Jacobian 
determinant  of  the  re-orientation  of  the  crystals  and  *A(s,0)  =  -4q(s)  is  the  ODF  associated  with  the 

reference  map  and  can  he  thought  of  as  the  initial  texturing  of  the  material.  The  ODF  is  the  microstructural 
feature  associated  with  every  material  point  in  a  simulation  of  industrial  forming  process.  The  two- length 
scale  sensitivity  framework  for  designing  deformation  processes  to  control  properties  is  shown  in  Fig.  7. 

The  mathematical  definition  of  the  optimization  problem  for  the  control  of  microstructure-sensitive  prop¬ 
erties  is  stated  as  follows: 


min r;  j8)  -  flrfos’red(i))2  dV  (3.2) 

where  Qdc3ired  js  the  desired  material  property  represented  at  a  material  point  on  the  final  product  and  (3 
is  the  design  parameter  to  be  selected.  The  optimization  problem  involves  the  gradient  of  the  property  fi 

with  respect  to  the  process  parameters,  represented  as:  VO  =  (^j-,  . gjjjj-).  The  needed  gradient 

information  for  optimization  analysis  will  be  computed  as: 

fl  2,.  ■  ■  ,A,,0,  ■  ■  ■  ,0,  A&.0, ...  ,0) 

A  & 


m 

dpi 


(3.3) 


Figure  7:  Description  of  the  parameter  sensitivity  analysis  involved  in  a  two- scale  deformation  problem.  The 
macro-process  design  variable  (s)  control  the  deformation  gradient  at  each  material  point  These  macro-fields 
define  the  underlying  microstruciure  and  properties  at  the  material  point  In  the  micro- sensitivity  analysis, 
we  compute  how  small  perturbations  in  the  deformation  gradient  at  the  material  point  affects  microstructure 
evolution .  On  the  macro- sensitivity  analysis ,  this  perturbation  to  the  deformation  gradient  at  each  material 
point  is  induced  by  the  perturbations  to  the  macro- design  variables.  Knowing  how  a  macro  desu/n  variable 
affects  microstructure- sensitive  properties  will  allow  us  to  provide  feasible  materials  design  obtainable  through 
processing . 

o 

where  ft  represents  the  sensitivity  of  ft  (with  respect  to  the  design  parameters  0).  Since  each  property  of 
the  poly  crystal  is  obtained  by  averaging  relations,  it  Is  essential  in  the  proposed  analysis  to  compute  the 
sensitivity  of  texture  evolution  with  respect  to  the  process  design  parameters  0. 

Computation  of  sensitivity  of  microstructure  would  involve  computation  of  sensitivity  of  evolution  of 
volume  fractions  of  different  crystal  orientation  components  over  the  deformation  history.  This  is  computed 
by  differentiating  the  DDF  evolution  equation  in  order  to  account  for  sensitivity  of  ODF  to  process  variables. 
Consider  the  Lagrangian  version  of  the  DDF  conservation  equation  (Eq,  (3.1)),  Design  differentiation  of  this 
equation,  assuming  that  the  initial  texture  is  independent  of  the  design  parameters  (/3),  leads  to  the  following: 


[  U(M;0,A/3)  =  -  f  MsJ\0)  J(sJ;0,A0)da  (3.4) 

Jn  J  n 

where  J  (sJ\0,&0)  —  det(Vr(s,t\ 0,  A/3))  —  J(s,t:0)  |v*  r  (s,  t:  0,  A/3)j .  Here,  sensitivities  of  the 

reorientation  vector  (r)  is  calculated  using  sensitivity  of  the  elastic  part  of  deformation  gradient  computed 
by  differentiating  parameters  in  the  constitutive  integration  algorithm  with  respect  to  process  parameters  of 
interest, 

Once  the  sensitivity  field  for  the  ODF  is  obtained,  it  is  straightforward  to  compute  the  sensitivities 
of  Cauchy  stress  for  use  in  the  macro-  sensitivity  equations  using  sensitivity-averaging  laws.  Evolution  of 
properties  such  as  yield  strength  and  Youngs  modulus  can  be  tracked  and  controlled  at  each  point,  in  the 
forged  component  by  knowing  the  sensitivities  of  these  parameters  to  the  process  variables. 

Described  below  is  the  analysis  for  the  development  of  a  total  Lagrangian  sensitivity  formulation  for  the 
kinematic  problem  at  the  macro-scale.  Let  the  reference  configuration  be  B®.  The  design  differentiation  of 


the  macro-scale  equilibrium  equation  results  in: 


V0.  { P )  +  f=  0  (3.5) 

where  (P)  is  the  poly  crystal  averaged  PK-I  stress.  A  variational  form  for  the  sensitivity  equilibrium  equation 

o  o 

(for  parameter  sensitivity)  can  be  posed  as  follows:  Evaluate  x=  x(Xyt;0tA0)  such  that 

f  (P)  *  VQ r}dVQ  =  f  \  .7}  dAo  (3.6) 

Bq  BBq 

for  every  a  kinematically  admissible  sensitivity  deformation  field  expressed  over  the  reference  configura¬ 
tion  E0  Details  on  the  methodologies  to  obtain  sensitivities  for  three-dimensional  die  and  preform  design 
problems,  and  to  model  sensitivities  of  die-workpiece  contact  through  design-differentiation  of  a  regularized 
contact  problem  can  be  found  in  our  work  in  [7].  Once  the  final  sensitivities  of  the  texture  to  process  variables 
is  found  at  the  end  of  simulation  of  a  deformation  process,  a  variety  of  gradient  optimization  algorithms  can 
be  Invoked  to  minimize  the  objective  function  defined  in  Eq.  (3.2)  for  tailoring  the  property  distribution  in 
the  material. 
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Figure  8:  Control  of  Youngs  Mod¬ 
ulus  variation  in  the  exit  cross- 
section  by  designing  optimal  die 
shapes,  (a)  shows  Youngs  Modulus 
distribution  at  the  end  of  the  ex¬ 
trusion  process  on  a  sample.  The 
polycnjstalline  texture  at  a  point  in 
the  exit  cjvss  section  is  indicated , 
(b)  Variation  in  the  Youngs  Modu¬ 
lus  at  the  exit  cross  section  for  dif¬ 
ferent  die  designs.  Optimal  design 
is  achieved  in  the  fourth  iteration. 


Here,  we  show  optimization  of  Youngs  Modulus  variation  in  the  exit  cross-section  during  metal  extrusion. 
Fig.  8  shows  Youngs  Modulus  distribution  at  the  end  of  the  extrusion  process  on  a  sample.  The  poly  crystal  line 
texture  at  all  points  in  the  exit  cross  section  is  controlled  by  modification  of  die  shapes.  The  efficiency  of  the 
methodology  is  exemplified  by  achievement  of  optimal  Youngs  Modulus  distribution  in  only  four  iterations. 
The  variation  in  the  Youngs  Modulus  at  the  exit  cross  section  for  different  die  designs  are  compared  in 
Fig.  8(b).  In  this  problem  the  objective  was  designed  so  as  to  minimize  the  variation  of  Youngs  Modulus 
over  the  exit  cross  section. 


4  Development  of  continuum  sensitivity  analysis  technique  for  de¬ 
sign  of  complex  three  dimensional  deformation  processes. 

We  have  extended  the  sensitivity  analysis  technique  to  allow  design  of  complex  3D  deformation  processes.  In 
extending  the  design  simulator  to  realistic  3D  simulations  some  important  issues  were  addressed.  For  complex 
dies,  as  is  usually  encountered  in  realistic  metal  forming  simulations,  die  surfaces  were  generated  using  finite 
element  meshes  followed  bv  the  construction  of  a  smooth  surface  from  the  individual  element  patches.  We 
employ  Gregory  patches  for  construction  of  a  smooth  surface  by  interpolating  the  finite  element,  nodes.  This 
also  provides  tangential  and  normal  continuity  between  adjacent,  patches.  In  this  work  automatic  hexahedral 


Figure  9:  Initial  preform  arid  final  forged  product,  for  the  steering  link  problem  (Example  6}. 

remeshing  was  implemented  using  the  meshing  software  CUBIT.  The  specific  remeshing  algorithm  that  was 
used  is  known  as  THEX.  In  this  approach  an  unstructured  hexahedral  mesh  is  generated  by  first  meshing 
the  workpiece  with  a  tetrahedral  mesh  following  by  division  of  each  tetrahedron  into  four  hexes.  To  aid  in 
speeding  up  the  solution  process  for  complex  forging  processes,  the  design  simulator  was  parallelized  using 
MPI.  The  simulator  was  dynamically  linked  to  the  parallel  toolbox  PetSc  for  parallel  assembly  and  solution 
of  linear  systems.  In  particular,  for  solution  of  linear  systems  a  GMRES  solver  along  with  block  Jacobi  and 
ILU  preconditioning  from  the  PetSc  toolbox  was  employed. 

Numerical  Example:  Isothermal  preform  design  for  a  steering  link 

The  example  deals  with  the  isothermal  preform  design  for  a  steering  link  to  ensure  the  required  shape  is 
obtained  while  the  material  wastage  is  minimized.  The  Howt  rule  for  the  alloy  is  given  by 

<«> 

where  eb  =  0.002  s~l,  s  —  150  MPa  and  n  =  5  while  the  elastic  parameters  are  taken  as  A  —  14423  MPa 
and  /z  —  9615.4  MPa. 

For  the  optimization  problem  we  start  with  an  initial  preform  shape  which  is  a  cylinder  of  much  smaller 
volume  than  the  volume  of  the  die  cavity  (Fig. 9).  The  preform  surface  is  discretized  using  two  design 
variables  which  are  essentially  the  radii  at  the  two  ends  of  the  cylinder  while  the  length  of  the  cylinder 
is  fixed.  The  forging  velocity  was  assumed  to  be  Q.Qimm/s  while  the  stroke  was  fixed  at  0,38  mm.  The 
total  simulation  time  was  38  seconds  and  partially  adaptive  remeshing  operations  were  performed  every  3 
seconds.  The  number  of  elements  varies  between  2CQ00  and  40000  during  the  remeshing  operations  within 
the  optimization  iterations.  Using  symmetry  only  one-eighth  of  the  actual  domain  evolution  is  simulated. 
The  BFGS  algorithm  was  used  for  the  gradient  optimization  process.  The  preform  shape  and  the  final  forged 
product  for  the  optimal  iteration  are  shown  in  Fig.  10. 

5  Development  of  stochastic  variational  multi-scale  model  with 
explicit  subgrid  modeling 

The  focus  of  this  work  is  to  model  uncertainty  induced  by  micro-scale  variability  on  the  macro-scale  outcomes. 
We  have  developed  mathematical  models  for  direct  incorporation  of  the  inherent  randomness  and  the  effect  of 
modeling  assumptions  in  the  design  of  upscaling  methods.  The  potential  applications  for  the  problem  include 
solving  PDEs  in  random  heterogeneous  microstructures,  which  is  a  key  problem  of  interest.  Heterogeneous 
microstructure  realizations  exhibit  property  variations  at  a  much  smaller  scale  compared  to  the  size  of 
the  computational  domain  D.  Performing  a  fully- resolved  calculation  on  these  mic restructures  becomes 
computationally  expensive.  We  consider  a  computational  scheme  that  involves  solving  for  a  coarse-soiution 
while  capturing  the  effects  of  the  fine  scale  on  the  coarse  solution. 


Figure  10:  Optimal  preform  and  final  forged  product  for  the  steering  link  problem  (Example  6) . 

Problem  definition  and  variational  formulation  Let  XL  T  and  0  denote  a  dosed  polygonal  domain,  a 
time  interval  and  a  suitable  probability  space,  respectively  The  transient  diffusion  in  a  heterogenous  medium 
with  a  spatially  varying  random  diffusion  coefficient  fc(x,u)  can  be  written  as  the  following  stochastic  partial 
differential  equation  (SPDE): 

utt  =  V  *  (jfc(x,  w)Vu)  +  / (xs  t,  w),  x  G  XL  t  G  T  and  bj  G  fi,  (5.1) 

where,  ult  denotes  the  partial  time  derivative  du/dt.  The  solution  u  =  u{x,f,u;)  and  the  source  term 
/ (x,  tyW)  are  real-valued  space-time  stochastic  fields.  For  closure*  we  also  assume  the  following  initial  and 
boundary  conditions: 

u{x,  £>u;)  —  w5(x,£,u>),  x  G  dV ,  t  G  T  and  (litQ,  (5,2) 

u(X]0ftu)  —  xGP,  and  uj  G  ft.  (5.3) 

The  variational  formulation  for  Eqs.  (5,1)— (5.3)  can  now  be  written  as:  Find  u  G  U  such  that  for  all 
v  G  V 

(u,t  ,v)  +  (fcVu,Vv)  =  (f,v),  (5.4) 

where 

( w,  u)  ^  f  f  uudxdXL  (5*5) 

Jn  Jt> 

Additive  scale  decomposition  and  variational  multiscale  method  In  the  variational  multiscale  ap¬ 
proach*  wre  consider  the  exact  solution  u  to  be  made  up  of  contributions  from  two  different  scales,  namely, 
the  coarse-scale  solution  uc ,  that  can  be  resolved  using  a  coarse- mesh  and  a  subgrid  solution  uF.  In  short, 
u  =  uc  +  uF .  The  main  idea  behind  the  VMS  (operator  upscaling)  method  is  to  develop  models  for  charac¬ 
terizing  the  effect  of  the  subgrid  solution  on  the  coarse-scale  solution  and  to  subsequently  derive  a  modified 
coarse-scale  formulation  that  only  involves  uc .  In  order  to  do  that*  we  split  the  variational  Formulation  given 
in  Eq.  (5  4)  as  the  following  coarse-scale  and  subgrid  equations:  Find  uc  G  Uc  and  uF  €  UF  such  that 

(■ uc,t  +uF,t ,  vc)  +  ( kVuc  +  WuF,  Vvc)  =  (/,  vc),  v/e  Vc,  (5.6) 

(uc>t+uF,t,vF)  +  (kVuc +  kVuFlVvF)  =  (f,vF),  V/e  VF.  (5.7) 

We  will  now  proceed  to  solve  Eq.  (5.7)  by  applying  localization  assumptions  to  obtain  an  approximate  model 
for  the  subgrid  solution  uF ,  Subsequently,  we  will  use  the  model  for  subgrid  solution  to  eliminate  in 
Eq.  (5.6)  to  obtain  a  modified  formulation  defined  only  in  terms  of  itc .  Subgrid  modeling  Assume  that 
the  spatial  domain  is  discretized  using  a  coarse- mesh  into  disjoint  sub-domains  T>^e\e  =  1,**,  ,Nd,  where, 
(e)  denotes  the  sub-domain  number.  We  will  refer  to  these  sub-domains  as  “coarse  elements”.  Let.  each 
coarse  element  be  further  discretized  using  a  subgrid  mesh  into  disjoint  sub-domains  (also  referred 

to  as  “subgrid  elements"). 


We  will  now  assume  that  the  subgrid  solution  is  a  sum  of  two  components  uF  anti  uF0  that  obey  the 
following  variational  equations: 

(uc,t)  vf  )  +  (uFlt  +  (W/,  VuF)  +  ( kVuF ,  Vt>F)  —  0}  (5,8) 

{ uF0,t,vF)  +  (kVuF0,VvF }  =  (f,vF),  (5.9) 

where,  uF  incorporates  entire  coarse-scale  solution  information  and  uF0  is  independent  of  the  coarse-scale 
solution.  The  dynamics  of  uF0  is  driven  by  the  projection  of  the  source  term  /(x,  £,uj)  onto  the  subgrid  trial 
function  space  VF . 

By  Eq.  (5.8),  the  term  iiF  behaves  as  a  mapping  from  coarse-scale  solution  to  the  subgrid  solution. 
Hence,  we  refer  to  uF  as  the  coarse- to-subgrid  (C2S)  map.  Owing  to  the  affine  nature  of  Eq.  (5.9),  we  shall 
refer  to  u* 0  as  the  affine  correction  term.  In  order  to  localize  the  calculations  of  the  subgrid  solution,  we 
will  consider  restrictions  of  Eqs.  (5.8)  and  (5.9)  to  each  coarse  element  K  The  derivations  in  subsequent 
discussion  are  performed  for  a  single  coarse  element 

C2S  map  and  niultiscale  basis  functions  Let  us  assume  that  in  a  fully- resolved  direct  numerical 
simulation*  the  dynamics  of  the  exact  solution  can  be  captured  using  a  fine  time-step  of  St.  Since  the  length 
scales  of  interest  in  the  coarse  solution  are  far  greater  than  the  smallest  length  scale  in  the  exact  solution,  we 
assume  that  the  coarse  time-step  A t  is  much  larger  in  comparison  to  St.  Let  us  consider  a  coarse  time-step 
At  —  Let  t*  be  the  local  time  coordinate  defined  such  that  at  tni  t*  =  0  and  at  tu+ 1,  f  =  At. 

Let  us  also  assume  a  piecewise  polynomial  finite  element  representation  for  the  coarse  solution  inside  a 
coarse  element  (please  note  that  (ej  is  suppressed  henceforth  to  simplify  notation): 

Nbf 

uc(x,  t,w)  =  u^(t',u))9p(x),  (5.10) 

3=i 

where,  Nbf  denotes  the  number  of  finite-element  shape  functions  (piecewise  polynomials)  defined  on  the 
coarse  element.  [8T  9] .  Let  us  further  assume  a  truncated  generalized  polynomial  chaos  expansion  with  Pc  T  1 
terms  for  each  of  the  coefficients 

Pc 

4.(*'>w)  =  c.(«).  (s.ii) 

Thus,  the  stochastic  finite- element  (spatial  finite- element  +  GPCE)  representation  of  the  coarse  solution  can 
be  written  using  local  time  coordinate  as 

Nbf  pg 

uc(x,t,w)  =  (5-12) 

0^1 3=0 


where,  denotes  the  nodal  solution  (each  node  has  Pc  +  I  degrees  of  freedom)  and  are  polyno¬ 

mials  from  the  hypergeometric  Askey  series.  The  form  of  Eq.  (5.12)  has  been  used  in  our  earlier  works  [8,  9], 
For  the  C2S  map  uF ,  we  seek  a  representation  similar  to  Eq.  (5T2); 

Nbf  Pc 

uF(x,t,u)  =  (5,13) 

0—\ 3=0 

The  fine-scale  variational  formulation  of  Eq.  (5,8)  can  be  re-written  after  substitution  of  Eqs.  (5.12)  and  (5.13) 
is  simplified  as  follows: 


(5.14) 


Without  loss  of  generality,  we  can  assume  the  following  representation  for  the  coarse-scale  nodal  solutions 
inside  the  coarse  time-step 

(5.15) 

where,  and  u^3  denote  the  nodal  coefficients  in  the  GPCE  of  the  coarse  solution  at  the  start  and  end 
of  the  coarse  time  step  {see  Fig.  11).  A(t ')  and  B(t')  are  special  positive  functions  that  obey  the  following 
relations: 

A(tf)  +  B(t')  =  1,  A( 0)  »  1,  A(At)  =  0,  B( 0)  «  0,  and  B{At)  =  1.  (5.16) 

Eq,  (5.15)  involves  representation  of  a  function  as  a  convex  combination  of  two  functions. 


1 


Figure  11:  A,  Schematic  of  the  time  integration  framework:  A t  is  the  coarse-time  step  and  t *  is  the  local 
time  coordinate.  The  integration  parameters  A(t')  and  B(t')  are  shown  in  the  figure.  Also,  and  are 
identified  as  the  coarse  solution  fields  at  the  start  and  end  of  the  coarse  time  step,  respectively.  B.  Schematic 
of  a  typical  coarse  element  sub-domain:  The  coordinates  normal  and  tangential  to  the  element  edges  are 
denoted  by  the  letters  n  and  7,  respectively. 


The  representation  given  in  Eq.  (5,15)  is  quite  general  and  incorporates  several  well-known  time  integra¬ 
tion  rules.  For  example,  «A(£#)  =  (At  —  f}/At  and  B(tf)  =  f /At  yields  a  backward- Euler  time  integration 
rule. 

From  Eqs.  (5, 13)  and  (5.15),  we  can  write  the  C2S  map  as  follows: 

Nbf  Pc 

u£s(x,i,(j)  =  [upsAit'^psix^t1  ,u)  +  u^aB(t')4>0a{x,t',w)^ .  (5.17) 

0=1 a=Q 


We  can  now  write  Eq.  (5.14)  as 


(5.18) 


Note  that  tlie  above  equation  is  fully  characterized  based  on  the  values  taken  by  Ups  and  and  the  subgrid 
basis  function  0^, 

We  are  looking  to  construct  a  localized  scheme  for  representation  of  the  subgrid  solutions.  Further,  this 
localized  scheme  should  hold  for  all  possible  values  of  the  coarse  nodal  coefficients  u^s  and  Hence, 
we  equate  the  terms  in  parentheses  to  zero  to  obtain  the  set  of  variational  formulations  defined  for  each 


combination  of  indices  0  and  s,  where,  0  =  l,,,.,NbF  and  s  =  Pc*  Now,  by  using  the  relations 

A(f')  +  B(t*}  =  1,  we  obtain  the  evolution  equation  for 

^(0^^  +  » v*  ^  +  <j>0a)t  VuF^  —  0.  (5.19) 

Affine  correction  term  The  affine  correction  term  uF0  as  defined  by  Eq.  (5,9)  leads  to  the  following 
strong  form  of  equations  inside  each  coarse  element  sub-domain 

uF0,t  +V  *  {kVuF0)  =  /,  x  €  D(c),  (5.20) 

Boundary  conditions  In  order  to  localize  the  computation  of  u  to  a  coarse  element,  we  need  an  ap¬ 
proximate  specification  of  its  boundary  conditions  along  coarse-element  edges.  A  schematic  showing  the 
calculation  of  boundary  subgrid  basis  functions  for  the  case  of  a  quadrilateral  coarse  element  is  shown  in 
Fig.  11.  In  the  VMS  simulation,  the  subgrid  basis  functions  are  generated  upto  a  cut-off  time  (referred  here 
as  the  bum-in  time)  that  is  chosen  such  that  the  effects  of  the  subgrid  component  of  the  initial  condition 
do  not  affect  the  coarse  solution.  This  also  ensures  that  the  subgrid  basis  functions  capture  sufficient  infor¬ 
mation  about  the  heterogeneities  at  the  subgrid  scale.  Note  that  the  affine  correction  term  does  not  have 
any  coarse-scale  solution  dependence,  hence,  the  assumption  that  this  term  goes  to  zero  on  coarse  element 
boundaries  is  justified  be,  uFD  =  0  on  coarse  element  boundaries.  Furthermore,  we  assume  that  the  affine 
correction  is  zero  at  the  start  of  each  coarse  time-step  he.  uFO(x,0,v)  =  0  (note  that  the  initial  conditions 
here  refer  to  the  local  time  coordinate  f'=0). 


Post- processing:  Fine-scale  solution  reconstruction  Let  us  assume  that  the  coarse-scale  solution, 
the  subgrid  basis  functions  and  the  affine- correct  ion  at  a  particular  time-step  are  uc(x,uj},  $0q(x,u>)  and 
aF0(x,u?),  respectively.  The  time  dependence  is  not  shown  here  for  clarity  sake.  Also,  we  emphasize  that  all 
calculations  in  this  section  are  performed  on  the  subgrid  mesh  associated  with  each  coarse  element  7>^e\ 
Inside  a  given  coarse  element  T> ^ ,  we  calculate  the  reconstructed  fine-scale  VMS  solution  as  follows; 

Nbf  Pc 

U  =  Y  E  +  hF0(x,w),  (5.21) 

P— 1 j=Q 

where,  is  the  sum  of  coarse  and  subgrid  basis  functions  defined  on  Let  us  now 

consider  the  GPCE  expansions  for  and  uF  0  as  follows: 

Pf  Pf 

=  Y  uFQ  =  Y  (5.22) 

r=Q  r=0 

wherein,  P?  +  1  is  the  number  of  terms  used  in  the  GPCE  expansion  of  and  uF0.  Now,  we  can  write 
Eq.  (5.21)  as  follows: 

Nbf  Pc  Pf  Pf 

U  =  £5I5Iu&®0Or(x)<r(w)<,(w)  +  U™(x)CrM ,  X  6  V{e) .  (5.23) 

0=1  3=0  r=0  r==0 


The  m-th  term  in  the  GPCE  of  the  reconstructed  fine-scale  VMS  solution  u  can  now  be  written  as  follows: 


um 


Nbf  Pc  Pf 


0=1 J=Q  r=0 


KrMGMCm  (<*>)) 

(C«(W)U(W)) 


+  «£>(*), 


(5.24) 


where,  (f(u>))  is  used  to  denote  the  mathematical  expectation  of  the  random  function  /(w). 

Numerical  Example:  Transient  diffusion  in  a  two- phase  microstructure  As  an  example,  a 
stochastic  variational  multi-scale  formulation  was  used  for  addressing  transient  diffusion  problems  in  random 
heterogeneous  microstructure  (  Fig,  12(a))  with  the  stochastic  parameter  being  the  diffusion  coefficient 


having  multiple  length  scales,  A  gray  scale  image  of  a  representative  two-phase  (a  and  /?)  microstructure. 
The  intensities  (I)  are  scaled  to  (0,  l]  ■  Let  pure  a -phase  and  pure  (3-  phase  be  associated  with  scaled  intensities 
7  =  0  and  7  =  1,  respectively. 

Probability  model  -  We  assume  that  the  pure  a  and  pure  ,3-phases  have  the  following  uniformly  distributed 
thermal  conductivities. 


fcnr(w)  =  ka o  +  kaiiliw),  and  fyj(w)  =  kpo  +  (5.25) 

where,  £i(u;)  and  are  two  independent  uniform  random  variables  defined  on  the  interval  [-1,1]*  We  use 
the  Following  mixture  model  for  defining  the  thermal  conductivity  at  a  given  location  on  the  microstructure: 

*(x>w)  =  (Mw)  -  ka(u))I(x)  +  ka(w),  (5.26) 

where,  /(x)  is  the  scaled  intensity  at  a  point  x  on  the  microstructure. 

GPCE  model  -  Since  the  input  distribution  is  uniform  (two-  dimensional),  a  two-dimensional,  third-order 
Legendre  chaos  expansion  (yielding  a  10  term  expansion)  was  used  for  representing  the  solution.  The  first 
few  terras  in  the  expansion  are  shown  below: 

«(M|0)  =  «0  +  «l€l(w)  +U2&(w)  +  «3(3?iM  -  l)/2, 

(W)6(u)  +  «5(3g(w)  -  l)/2  +  . . . ,  (5.27) 

The  main  idea  here  was  to  calculate  the  diffusion  dynamics  on  a  highly  coarse-mesh  by  including  localized 
solutions  to  sub-grid  problems.  The  method  captures  the  mean  field  (Fig,  12(b))  as  well  as  statistical 
quantities  such  as  the  first  order  fluctuations  (Fig.  12(c))  in  the  temperature  field  in  a  single  simulation. 


Figure  12:  Modeling  diffusion  in  heterogeneous  random  microstructures  using  a  variational  multi-scale  al¬ 
gorithm:  (a)  Microstructure  image  realization,  (b)  Mean  temperature  field  (coefficients  uQ  of  GPCE),  (c) 
First-order  statistics  (coefficient  ul  in  the  GPCE). 


6  Development  of  maximum  entropy  techniques  for  modeling  topo¬ 
logical  uncertainties  in  polycrystalline  metallic  microstructures 

Probability  distribution  functions  (PDFs)  providing  a  complete  representation  of  microstructural  variability 
in  3D  polycrystalline  materials  using  limited  information  is  difficult  to  obtain  since  this  inverse  problem  is 
highly  ill-posed.  We  use  the  maximum  entropy  (MaxEnt)  principle  to  compute  a  PDF  of  microstructures 
based  on  given  information  about  a  microstructural  system.  Microstructural  features  are  incorporated  into 
the  maximum  entropy  framework  using  data  obtained  from  experiments  or  simulations.  Statistics  of  material 
properties  are  then  extracted  by  interrogating  sample  microstructures  obtained  from  the  resultant  probability 
distribution  function  (PDF).  These  property  statistics  are  important  because  they  give  a  numerical  estimate 


of  how  the  material  may  perform  when  subjected  to  extremal  conditions,  A  particular  problem  of  interest 
addressed  here  is  the  determination  of  effective  behavior  of  three-dimensional  (3D)  poly  crystal  line  materials 
based  on  uncertainties  induced  due  to  randomness  in  grain  distribution  and  texture  of  microstructures. 

Some  important  concerns  which  are  addressed  in  this  work  are:  (i)  how  to  quantify  uncertainties  in  the 
micro-scale  based  on  limited  observable  behavior  at  the  macroscale,  (ii)  how  to  model  poly  crystal  microst  nur¬ 
tures  and  obtain  microstructural  samples  incorporating  the  uncertainties  computed  above,  and  (iii)  how  to 
obtain  statistics  of  homogenized  plastic  properties  from  this  statistical  description  of  microstructures.  In 
obtaining  the  distribution  from  which  microstructures  are  sampled,  we  use  the  principle  of  maximum  en¬ 
tropy  (MaxEnt).  MaxEnt  Is  a  popular  technique  that  is  utilized  to  construct  the  PDF  for  stochastic  fields 
when  limited  information  is  available  about  them.  Availability  of  limited  information  is  often  the  case  due 
to  high-costs  involved  in  performing  a  large  number  of  experiments  or  simulations.  For  instance,  grain  size 
distributions  of  polycrystalline  materials  cannot  be  experimentally  measured  easily  since  they  vary  even 
at  a  material  point  for  different  samples  of  the  material.  Such  serious  limitations  on  the  measurability  of 
stochastic  fields  has  made  MaxEnt  an  immensely  popular  choice  amongst  statisticians.  To  make  this  method 
feasible,  information-theorists  developed  a  mathematical  tool  for  measuring  entropy  of  systems.  Given  a 
random  variable  x  which  takes  discrete  values  given  by  xi,x2t . .  ♦ ,  xn ,  the  informational  entropy  of  x  is  given 
by 

n 

H  (p)  =  -^p(Xj)ios(p{Xi))  (6.1) 

i=  1 

where  p[Xi)  denotes  probability  of  the  discrete  event  xt  occurring.  Suppose  we  have  insufficient  knowledge 
about  the  distribution  p.  MaxEnt  provides  a  rationale  to  obtain  the  entire  probabilistic  variability  about 
the  variable  x.  It  is  to  be  noted  that  the  entropy  function  H(p)  is  convex  and  in  an  unconstrained  problem, 
it  achieves  its  maximum  value  when  all  the  possible  events  x*  are  equiprobable.  This  occurs  when  all  of 
the  p(x*}  are  equal.  We  call  this  particular  choice  of  p(x{)  as  g(x^),  The  distribution  g(xj)  does  not  lead  to 
maximum  entropy  when  there  is  some  knowledge  about  the  system  x,  which  denotes  the  vector  of  variables 
Xj.  To  distinguish  a  random  variable  from  a  random  field,  we  use  the  notation  x  for  the  latter.  Typical 
knowledge  about  microstructural  systems  of  polycrystalline  materials  is  given  in  the  form  of  lower-order 
statistical  moments  of  grain  sizes  as  well  a_s  average  texture  values  in  the  form  of  an  orientation  distribution 
functions  (ODF).  These  known  moments  are  denoted  as  A/*.  Features  such  as  grain  sizes  and  ODFs  are 
denoted  as  /*.  Knowledge  about  microstructural  features,  which  can  be  computed  using  experiments  or 
simulations,  is  posed  as  a  constraint  for  the  MaxEnt  optimization  problem  which  is  mathematically  stated 
as  follows: 

Find  :  p*(x)  —  argmaxH(p) 

p 

with  the  given  constraints 


Mi 

=  E(/,(x)) 

M2 

“  E(/3(x)) 

Mn 

=  E(/n{x)) 

(6.2) 

where  E  denotes  expectations  over  the  space  of  random  variables.  When  considering  topological  uncer¬ 
tainties  within  a  microstructure  (Problem  1  in  Section  5),  we  define  the  operator  E  for  grain  size  feature 
over  the  set  of  grains  within  a  microstructure  as  given  below: 

1  G 

=  (63) 

where  G  is  the  number  of  grains  in  the  microstructure  and  /  represents  the  grain  size  feature.  The  tUl 
moment  of  a  specific  feature  ft  is  computed  as  E(//).  The  distribution  that  is  finally  computed,  p*(x) 


is  the  most  uniform  distribution  satisfying  the  given  constraints  in  Eq,  (6.2).  This  means  that  p*  (x)  is 
closest  in  distance  to  q (x)  while  satisfying  the  available  information  about  the  system  given  by  Eq,  (6.2).  A 
popular  choice  of  distance  measure  is  the  Kullback-Leibler  (KL)  divergence  measure.  For  any  two  arbitrary 
probability  distributions,  m  and  n,  the  KL  divergence  is  given  as  follows: 

I<L(m\\n)  =  (6-4) 

Hence,  the  distribution  (x)  satisfies  the  additional  property  that: 

p* (x)  =  arg min  /iL(p||g)  (6.5) 

p 

while  satisfying  Eq.  (6.2).  This  KL  divergence  measure  is  crucial  in  reconstructing  microstructures  based  on 
grain  size  distributions.  The  problem  of  computing  a  probability  distribution  p(x)  given  the  constraints  in 
Eq.  (6.2)  is  ill-posed.  MaxEnt  is  a  technique  that  is  able  to  establish  a  unique  PDF  p*(x)  that  satisfies  the 
given  information.  The  maximum  entropy  technique  makes  the  problem  well-posed  by  imposing  additional 
restrictions  on  the  computed  PDF,  Utilizing  the  microstruct ural  features  such  as  grain  sizes  and  ODFs  (/,), 
the  PDF  that  is  obtained  using  the  M axEnt  technique  is  of  the  following  form: 

g  “  E  n-m  l  /«  (X) 

where  Xt  represent  the  Lagrange  multipliers  that  account  for  the  N  system  constraints.  This  exponential 
form  as  seen  in  Eq.  (6,6)  is  a  very  central  concept  in  the  MaxEnt  formalism,  A  number  of  numerical  options 
are  available  to  compute  a  microstructural  system  that  has  maximum  entropy  whthin  its  features  while 
simultaneously  satisfying  given  information  about  the  system.  A  Gibbs  sampler  is  used  that  starts  from  a 
random  micros  true  ture  and  after  a  period  known  as  burn-in,  generates  microstructures  that  are  sampled 
from  PDFs  of  the  form  given  in  Eq,  (6.6),  Samples  generated  after  the  burn-in  period  are  accepted  as 
microstructural  samples  from  the  desired  distribution.  The  algorithm  for  grain  sizes  follows  closely  the 
algorithm  given  in  [10,  11]. 

Numerical  Example:  Variability  in  Aluminum  polycrystals  driven  by  mean  grain  sizes: 
Literature  is  replete  with  many  instances  of  phenomenon  such  as  Hall- Fetch  effect  that  depend  only  on  the 
mean  grain  size.  In  such  a  situation,  the  variability  of  grain  sizes  within  a  microstructure  becomes  superficial 
since  clearly,  different  distributions  may  represent  the  same  mean  grain  size.  Hence,  we  extract  a  different 
feature  which  quantifies  how  mean  grain  sizes  vary  across  different  samples.  The  MaxEnt  distribution  is 
defined  over  a  set  of  microstructures,  each  characterized  by  its  mean  PDF,  Potential  applications  of  this 
example  relates  to  predicting  statistical  behavior  of  random  polycrystals  where  the  behavior  is  a  function  of 
mean  grain  sizes. 

This  problem  is  driven  by  variability  in  mean  grain  sizes  of  different  given  Aluminum  microstructural 
samples.  The  variability  in  mean  grain  sizes  as  well  as  textures  are  utilized  to  define  a  PDF  of  microstructures. 
Based  on  samples  from  this  PDF  and  using  a  homogenization  scheme  [6],  we  compute  bounds  in  the  plastic 
stress  strain  curve  of  Aluminum.  Using  input  microstructures  generated  (Fig.  13)  from  3D- phase  field 
simulations,  we  compute  the  expected  grain  sizes  (volumes)  of  each  microstructure. 

Lower  moments  until  the  fourth-order  of  these  (mean)  grain  sizes  are  extracted  and  a  maximum  entropic 
distribution  is  computed.  We  use  just  five  input  microstruct  ures  because  in  practice,  it  may  prove  costly 
to  compute  experimental  microstruct  ures  at  different  material  points.  Hence,  t  lie  input  information  relies 
only  on  the  limited  set  of  microstructural  samples.  The  resultant  distribution  based  on  two,  three  and  four 
mean-grain  size  moments  is  depicted  in  Fig,  14.  Each  point  in  the  ;r  —  axis  corresponds  to  a  microstructure 
with  the  corresponding  value  as  the  mean  grain  size.  Hence,  the  PDF  just  represents  the  probabilities  of  the 
cor  res  p  ond  i  ng  microst  r  ue  t  u  res , 

The  MaxEnt  distribution  was  computed  using  the  stabilized  version  of  the  conjugate  gradient  algorithm. 
The  computed  distribution  also  shows  that  the  randomness  induced  due  to  MaxEnt  reduces  as  we  incorpo¬ 
rate  more  information  about  the  source  of  randomness  (incorporated  using  higher  moments).  We  compute 
voronoi-cell  microstruct  ures  which  are  sampled  from  the  MaxEnt  distribution. 


p'm=j 


Figure  13:  Two  out  of  the  five  samples  of  microstructures  computed  using  phase  field  technique  that  was 
used  in  problem  L 


Figure  14:  A  MaxEnt  distribution  of  mean  grain  sizes  based  on  inputs  from  3D  phase-field  simulations. 


Further,  samples  from  the  tail  of  the  MaxEnt  distribution  were  computed  to  obtain  bounds  on  the  plastic 
properties.  Samples  near  the  tails  are  computed  as  follows:  Start  randomly  sampling  from  the  MaxEnt 
distribution  till  a  sample  which  falls  outside  the  05%  confidence  region  is  reached.  95%  confidence  region  is 
the  region  where  95%  of  random  samples  are  expected  to  lie  and  this  is  distributed  evenly  on  either  side  of 
the  mean.  This  is  accepted  as  a  tail  sample.  We  choose  one  tail  sample  from  each  extreme  of  the  MaxEnt 
distribution.  Each  sample  computed  from  the  PDF  corresponds  to  a  microstructure  with  the  corresponding 
mean  grain  size.  Some  microstructures  that  were  sampled  from  the  tail  are  shown  in  Fig.  15,  These  were 
meshed  using  the  commercially  available  software  CUBIT  using  hexahedral  meshes. 

Filially,  a  rate- independent  homogenization  method  is  utilized  to  compute  bounds  of  the  non-linear  plastic 
properties  of  the  reconstructed  microstructures  based  on  a  simple  tension  test  using  a  strain  rate  of  0.01  per 
second.  Bounds  in  the  properties  based  on  a  set  of  25  samples  are  shown  in  Fig.  1G> 

The  variability  results  due  to  difference  between  mean  grain  sizes  of  these  microstructures.  Tails  of  this 
distribution  correspond  to  extremes  of  mean  grain  size  and  were  crucial  in  defining  bounds  on  the  material 
behavior  of  Aluminum  polycrystals.  We  also  sawr  how  the  incorporation  of  different  amounts  of  information 
resulted  in  stricter  bounds.  Mathematical  characterization  of  such  distributions  between  microstructures 
pose  potential  applications  in  stochastic  simulations  such  as  SSFEM  which  are  still  driven  by  analytical 
assumptions.  Further,  statistical  characterizations  of  phenomenon  such  as  Hall-Petch  effect  can  be  made 
using  the  variability  in  mean  grain  sizes. 


Figure  15:  Samples  obtained  at  locations  A  and  B  of  Fig.  14. 


Figure  16:  Bounds  of  the  material  strength  computed  from  microstructures  sampled  near  the  tail  of  MaxEnt 
distribution. 
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New  Discoveries 

(a)  The  first  stochastic  variational  multi-scale  algorithm  with  explicit  subgrid  modeling  reported;  (b)  The 

first  techniques  reported  For  full- probabilistic  modeling  and  design  of  deformation  processes  (c)  Multi-scale 
computational  simulator  for  design  of  mechanical  properties  such  as  stiffness  and  strength  during  deformation 
processes  by  tailoring  underlying  microstructures  (d)  New  techniques  for  reconstruction  of  3D  microstructures 
from  limited  statistical  information  (e)  Techniques  for  computing  property  variabilities  due  to  topological 
uncertainties  in  the  microstructure. 


